Coarse-grained loop algorithms for Monte Carlo simulation of quantum spin systems 
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Recently, Syljuasen and Sandvik proposed a new framework for constructing algorithms of quan- 
tum Monte Carlo simulation. While it includes new classes of powerful algorithms, it is not straight- 
forward to find an efficient algorithm for a given model. Based on their framework, we propose an 
algorithm that is a natural extension of the conventional loop algorithm with the split-spin repre- 
sentation. A complete table of the vertex density and the worm-scattering probability is presented 
for the general XXZ model of an arbitrary S with a uniform magnetic field. 
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I. INTRODUCTION 

Among many numerical techniques for condensed mat- 
ter physics, the Monte Carlo method is a popular choice 
when a long correlation length or a small excitation gap 
is anticipated. Apart from the negative sign difficulty, 
most of the shortcomings of the quantum Monte Carlo 
(QMC) method for spin systems have been removed or 
reduced. In particular, the QMC for finite temperature 
based on the path-integral representation has been im- 
proved considerably during the past decade. The im- 
provement was achieved mainly by the development of 
the loop-cluster algorithms jl) and related methods. For 
instance, the critical slowing down was tamed by the 
loop-cluster algorithms for a broad class of quantum spin 
systemsj2]. It was shownjH that the typical size of the 
clusters coincides with the correlation length. Because 
of this property, an effective update of configurations is 
possible. The other slowing down, due to small inter- 
vals for discretization of the imaginary time, was com- 
pletely removed also by the loop-cluster QMC[Q, which 
became even more evident by the extension to continu- 
ous imaginary time[^. An efficient measurement of some 
of important off-diagonal quantities was made possible 
through the improved estimator 

One of the difhculties that have been left unsolved 
until recently was the freezing of configurations due to 
an external field competing with the exchange couplings. 
In the conventional framework of the loop-cluster algo- 
rithms, the field term does not affect the graph assign- 
ment probabilities. It is taken into account only in the 
flipping probabilities of clusters. As a result, the clus- 
ter size does not correspond to the physical correlation 
length any more. It was demonstrated that this diffi- 
culty can be removed, in the case of the 5 = 1/2 antifer- 
romagnetic Heisenberg model, by introducing two singu- 
lar points at which the local conservation rule of particle 
number (or magnetization) is violated. These singular 
points are called "worms." This extension of the config- 
uration space makes it possible to take the external field 



into account in the hopping probability of worms. 

Another difficulty is large memory requirement due 
to the split-spin representation^. When one uses the 
loop algorithm for a spin problem with large S, it is cus- 
tomary to replace each spin operator by a sum of 25* 
Pauli matrices. Therefore, for larger S, the algorithm 
consumes more memory. The stochastic series expansion 
(SSE)||, ^ does not have this difficulty, since it works di- 
rectly on the original spin-configuration space. The SSE 
is based on the high-temperature series expansion of the 
partition function, rather than the path-integral formula- 
tion. However, it was pointed outjl^, |l^ that these two 
apparently different formulations are essentially equiv- 
alent in the limit of the infinite order expansion. The 
apparent difference was due to the different updating 
method, rather than the formulations themselves. 

Quite recently [Tc|], Syljuasen and Sandvik introduced 
the notion of "directed loops" and proposed a framework 
that accommodates all of the above-mentioned ideas, i.e., 
the loop updating, the worm updating, and the two for- 
mulations. Their framework can be compared with Kan- 
del and Domany's framework [|3) for the loop-cluster al- 
gorithms. In fact, the mathematical formulation of the 
Syljuasen-Sandvik (SS) scheme has a similar structure to 
Kandel and Domany's (see Appendix^, and the result- 
ing algorithm coincides with the loop-cluster algorithm 
in some cases. In this sense, the SS scheme can be viewed 
as a generalization of the Kandel-Domany framework. 

An algorithm based on the SS scheme is characterized 
by the scattering probabilities of worms. Although the 
detailed balance condition imposes a set of equations to 
be satisfied by these probabilities, there are still a lot 
of degrees of freedom. Some of the solutions to these de- 
tailed balance equations lead to the single-cluster version 
of the conventional loop-cluster algorithm at zero mag- 
netic field, which are known to be efficient. However, 
all the solutions are not necessarily efficient or practi- 
cal. There are obviously many bad solutions in which the 
"back-tracking" probability |lO| or "turning-back" prob- 
ability are dominating. In addition, the straightforward 



solutions of the heat-bath type do not work either, as we 
see below in Sec. 0. Although an efficient solution was 
discussed [p^ for S' = 1/2 XX Z models, the prescription 
was not given for general S. 

Similar to the Kandel-Domany framework, the SS 
scheme does not give a concrete prescription for obtain- 
ing a good solution that leads to an efhcicnt algorithm for 
specific models. A rule of thumb for obtaining a good so- 
lution is to minimize the turning-back probability. How- 
ever, even if the turning-back probability is fixed, we still 
have many degrees of freedom to play with, and the ef- 
ficiency of the algorithm strongly depends on the choice 
of the worm-scattering probabilities, as we demonstrate 
in Sec. 0. While this freedom can be quite useful for 
constructing new types of efficient algorithms, it makes 
finding a reasonable solution a nontrivial task. 

In this paper, we propose a natural extension of exist- 
ing algorithms that determines a unique set of scattering 
probabilities of worms. The resulting algorithm is within 
the SS scheme and expected to be efficient for a wide 
class of quantum spin systems. The algorithm can be 
obtained by the coarse graining mapping applied to an 
algorithm in the split-spin representation. In Sec. |l[ we 
first discuss algorithms in the split-spin representation. 
In Sec. [II, we show how we can obtain an algorithm 
in the original spin representation by coarse-graining the 
split-spin algorithm. Based on this general prescription 
for the coarse-grained algorithms, we present in Sec. |^ 
a complete table of the worm scattering probability for 
the XXZ models with arbitrary S. Finally, in Sec. 
we compare the present algorithm with other algorithms 
such as the directed loop algorithm with the straightfor- 
ward heat-bath solution. 



II. ALGORITHMS IN THE SPLIT-SPIN 
REPRESENTATION 

It has been pointed out in previous papers (l^, ^ that, 
in the limit of the infinite order expansion, a Monte Carlo 
algorithm based on the series expansion can be reformu- 
lated in terms of the language of the path-integral with 
continuous imaginary time, and vice versa. In what fol- 
lows, therefore, we describe algorithms in this limit for 
the sake of simplicity, and use the path-integral language. 
The translation into the series-expansion language and 
its modification for a finite order expansion should be 
straightforward. 

A simulation based on the path-integral representation 
(or the SSE in the infinite order expansion limit) can be 
visualized in a (d + 1) - dimensional space-time where d 
is the real-space dimension. At each point in this space- 
time an integral variable is defined, and it takes on one 
of the 2S + 1 values —S, ~S -t- 1, . . . , and S. In the case 
of S' = 1/2 the variables are one-bit (or Ising) variables. 
Accordingly, we consider world lines, which are trajec- 
tories of up-spins in this space-time. In the present pa- 
per, we use the term "world-line configurations" to refer 
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FIG. 1: Various objects in (1 -|- l)-dimensionaI space-time 
in the case of S = 1. The vertical direction corresponds to 
the temporal coordinate whereas the horizontal direction cor- 
responds to the spatial coordinate. Kinks {Ki,K2, and K3), 
vertices (¥1,^2, V3, and V4), and worms (Wi and W2) are 
shown. A number (0, 1, or 2) printed on every segment means 
the number of particles [= (r) -I- S]. Horizontal lines rep- 
resent vertices. The kinks coincide with the vertices at which 
particles jump from segments to segments. 



to the spin configurations in the space-time for general 
S, although they are not represented by simple lines for 
S" > 1/2. A Monte Carlo algorithm is nothing but a pro- 
cedure by which the world-line configuration is updated 
so that the limiting probability distribution may coincide 
with the weight of the configuration, i.e., the exponential 
of the action. 

In the SS scheme, we deal with objects defined in the 
(d + 1)- dimensional space- time (Fig.|l|). A vertical line 
of length f3 represents a spin. A kink is a point at which 
the local spin configuration changes. A particle (or an 
up-spin) jumps from one vertical line to another only at 
kinks. In models without particle number (or magneti- 
zation) conservation, a point at which a particle disap- 
pears or appears is also a kink. Every kink is located 
on a vertex. Vertices in the SS scheme play a role com- 
parable to that of local graph elements in the conven- 
tional loop-cluster algorithms. In particular, for models 
in which the magnetization conserves, vertices are rep- 
resented by short horizontal lines, each connecting two 
or more neighboring vertical lines. If a vertex connects 
two lines, we call such a vertex four legged since it joins 
four segments, where a segment is a part of a vertical line 
which is delimited by two vertices. A worm is a kink of 
a special kind located on a segment jll]. A worm can 
move continuously as the simulation proceeds, while lo- 
cations of ordinary kinks and vertices are fixed until they 
are deleted. In the applications discussed in the present 
paper, there are only zero or two worms at the same time 
in the whole system. 

For quantum spin systems, one cycle of update in the 
SS scheme consists of the following operations on these 



objects: 1) assigning vertices to a given world-line con- 
figuration, 2) creating a pair of worms, 3) letting one of 
them move along segments and be scattered by vertices 
until it comes back to the other worm to be annihilated, 
and 4) deleting all the vertices with no kinks on them. 
In the rest of the present paper, we see these operations 
in more detail. 

When a world-line configuration is given, we first as- 
sign vertices. Vertices are assigned to every part of the 
system probabilistically with a density that depends on 
the local world-line configuration. In addition, all the 
kinks are regarded as vertices. After placing all vertices, 
we choose a point on a segment at random and create 
a pair of worms there. Then, one of the worms starts 
moving. As a worm passes a point in the space-time, 
it changes the local spin value there. When the moving 
worm encounters a vertex, it may be scattered. The out- 
going direction after the scattering is determined stochas- 
tically with certain predetermined scattering probabili- 
ties. When a moving worm meets the other worm, they 
annihilate. Therefore, what we have to specify in order 
to define an algorithm are the density of vertices and the 
scattering probability of worms. The SS scheme imposes 
conditions on these. The conditions are summarized in 
Appendix 

When spins in a given model are larger than S=l/2, 
it is customary to replace the original spin operators by 
the sum of 25* Pauh spins 0, i.e.. 



(a) 



(b) 



Qa nc 



2S 



where af^ is an S 



1/2 spin operator. The partition 
function is expressed in terms of these a degrees of free- 
dom. Since the original phase space corresponds to the 
subspace in which 

holds, we have to project out all the states orthogonal 
to this subspace to obtain the correct partition function. 
This can be done by inserting the projection operator P: 



Z = Tr{,.j(e-'^^[{^'>l) = Tr{..,.j(Pe 



-/3H[{S.}] 



Hereafter, the representation based on Si degrees of free- 
dom is referred to as the "original spin representation" 
whereas that based on (7^^ the "split-spin representa- 
tion." 

For many models, it is rather straightforward to ob- 
tain an algorithm in the split-spin representation. For 
example, we can obtain an algorithm for the XX Z mod- 
els with arbitrary magnitude of spins S from that for the 
corresponding model with S = 1/2, simply by regarding 
the former as a superposition of many of the latter. If 
we do so, we consider 2S vertical lines for each original 
spin. Accordingly, a point in the space-time is specified 
by three numbers ((z,^),t) rather than two (?,t). The 



FIG. 2: Two types of 12-legged vertices that appear in the 
S — 3/2 SU(4) models in the split-spin representation. 



coupling between two original spins Si and Sj is trans- 
formed into (25*)^ couplings, each couples CTj.p and cr^ 

-H,, = JS^S^ + JSfS^ + J'S!S] + ^{S:^ + S]), 



where J > 0, h = Hp/{2S), and Hp is the external field 
per bond (e.g., Hp — H/d for the hypercubic lattice 
where H is the external field per site). For the vertex 
assignment, we apply the procedure for the directed loop 
algorithm for S = 1/2|10 to every one of (25')^ combina- 
tions of a spins. To be more precise, the density for the 
vertex between two split spins is the same as that in the 
directed loop algorithm for S* = 1/2 with Hp replaced by 
h = Hp/{2S). Similarly, the worm-scattering probabili- 
ties for S" = 1/2 can be used for split spins with the same 
modification of Hp. 

For the projection operator, we do essentially the same 
as we usually do in the conventional loop algorithm for 
S > 1/2 iQ. In the present framework, we represent it 
by special vertices, each located at r = /? connecting all 
the {2S) vertical lines on a site i. To be specific, when 
a worm moves upwards along the vertical line (z, /i) and 
hits the point ((i, n), /3) from below, it jumps to ((i, v),0) 
and go on upwards. The line to which the worm jumps, 
is chosen with equal probability among those on 
which the local spin state is the same as the spin state 
right above the incoming worm. Namely, it is chosen 
among such v's that ai^{0) 



(Tif^iP) may hold. 



III. COARSE GRAINING 

One of the drawbacks of the split-spin representation 
mentioned above is that it may require much more mem- 
ory than the original spin representation. For example, 
in the loop algorithm for the SU {N) models [|l5| , we in- 
sert graphs that involve all Pauli spins on two neigh- 
boring sites at the same imaginary time (Fig.|^) . In the 
split-spin representation, insertion of a graph of this type 
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creates 2(iV — 1) = 45* new segments. Since the mem- 
ory requirement is roughly proportional to the number 
of segments, a loop algorithm for the SU (N) model re- 
quires memory resources proportional to 2(iV — 1). If we 
can construct an algorithm in the original spin represen- 
tation, insertion of a graph would create only two new 
segments. This leads to a memory requirement smaller 
by factor 1/{N — 1) than that of the split-spin represen- 
tation. 

Another drawback is the lack of portability of the code. 
In the split-spin representation, there are many kinds of 
vertices, in principle, depending on the number of legs. 
Therefore, we have to change the core part of the code to 
accommodate new kinds of graphs for each model unless 
we implement all possible sorts of graphs initially, which 
is impractical. On the other hand, in the coarse-grained 
representation, all the vertices are four-legged (for mod- 
els with two-body interactions) and there are only four 
different types of scattering of worms. Therefore, the 
core part of codes for all models are the same except for 
the densities of vertices and the scattering probabilities 
of worms. For example, if we have a code based on the SS 
scheme for the SU{N) model we can immediately obtain 
a code for the XY model simply by changing the arrays 
of the probability tables. 

In order to take full advantage of the SS scheme, there- 
fore, we have to construct probability tables for algo- 
rithms based on the original spin representation rather 
than the split-spin representation. For this purpose, we 
consider a "coarse-graining" map and its stochastic in- 
verse. The map is basically disregarding the detailed 
information of split spins. The inverse of the map is to 
choose stochastically one of the split-spin configurations 
which are transformed by coarse graining into a given 
original-spin configuration. With these maps and an al- 
gorithm in the split-spin representation, we can construct 
an algorithm in which we manipulate only original spin 
degrees of freedom. 

To illustrate the idea in more detail, we again take 
the XXZ model with an arbitrary magnitude of spins 
S. We can define a coarse-graining map from a split-spin 
world- line configuration S into an original spin world-line 
configuration 5' as 

S = {ai^{T)}^S = {Si{T)}, 

where (Jinir) is the value of af^ at the imaginary time 
T, whereas Si{T) is the sum of them, i.e., Si{T) = 
^^(Tip(r). Similarly, we can define a map for vertices. 
Since the only interaction is of second order in the spin 
operators, a vertex of the XXZ model has only four legs 
(i.e., connects only two lines). If a vertex connects two 
lines (i,/i) and we associate with it a vertex that 

connects two coarse-grained lines i and j. Of course, in 
the latter representation, the information concerning the 
indices /z and v is missing. 

Obviously, these mappings are many-to-one mappings. 
However, wc can define the inverse of this coarse-graining 
map. In the following, we adopt the "particle" picture 



in which an up-spin is regarded as a particle whereas 
a down-spin a hole. Correspondingly, wc use particle 
numbers li{T) = 0, 1, . . . , 25* and nifj_{T) = 0, 1, instead of 
Si{T) and (Jifi{T), to specify local states of spins. These 
are related to S'j(r) and crj^(T) by Si{T) = li{T) — S and 
(^inir) = niu,{T) - 1/2. 

The inverse mapping of a local state is rather simple. 
Suppose that a model is an = 1 model and a local spin 
state at the point of interest is Z,(t) = 1 in the coarse- 
grained representation. There are two split-spin states 
that are mapped to this state, i.e., {nii{T),ni2{T)) = 
(1,0) and (0,1). Both configurations are chosen with 
the same probability (i.e., 1/2) since there is no reason 
to put any bias. For general S, all configurations that 
satisfy = / are chosen with the same probability, 

where / = 0, 1, . . . , 25 is the local state on the coarse- 
grained line. 

The inverse mapping of a vertex can be defined in a 
similar way. When two space-time points {i, r) and {j, r) 
are connected by a vertex (with no kink on it) in the 
coarse-grained representation, we can map it to a vertex 
connecting ((i,/Lt),T) and {{j,v),T) with some probabil- 
ity. When S = 1, there are four different ways of choosing 
H and v that arc to be connected. However, the proba- 
bility for taking one of them is not 1 /4 in this case. This 
is because the density of vertices depends on the spin 
states at their legs. If, for example, the density of ver- 
tex between two particles is higher than that between a 
particle and a hole, a given coarse-grained vertex should 
be mapped to the former with larger probability than 
the latter. In other words, the probability for associat- 
ing a coarse-grained vertex with a particular split-spin 
vertex should be proportional to the density of the lat- 
ter. Therefore, the probability for associating a coarse- 
grained vertex between i and j with a split-spin vertex 
between (i,/x) and (j, ly) is given by 



Here, p''^^, is the density of vertices in the split-spin rep- 
resentation where the local spin values are n and n' at 
the legs of vertices. 

The coarse-graining map and its inverse can be used 
for obtaining the vertex density and the worm-scattering 
probability in the coarse-grained representation from 
those in the split-spin representation. Suppose an 
imaginary-time interval in which the state of two neigh- 
boring sites i and j arc specified by / and rrj, respectively. 
The site i consists of I particles and I (= 25* — 1) holes 
whereas the site j consists of m particles and rh holes. 
Then, there are Zm, Ifh, Im, and Ifh possible combinations 
of 11, 10, 01, and 00 pairs of a spins, respectively. Since 
we assign a vertex with density p^^^ for each 11 pair, the 
total density of vertices connecting 11 pairs is Imp'ff'. 
The densities for other combinations can be obtained in 
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FIG. 3: Scattering at a vertex for S = 3/2 in the split-spin 
representation. The thick horizontal lines represent vertices. 
The spin-lowering worm is indicated by an open triangle. A 
solid line represents a world line with ni^ij) — 1 whereas 
a dotted one represents ni^ij) — 0. In case (a), the worm 
cannot be scattered. There are two cases [(b) and (c)] where 
the worm can be scattered. 



a similar fashion. Thus, the total density of vertices is 
pim = Impfl^ + lfhp''^Q^ + Imp'^fl^ + Ifhp'^fgK (2) 



for two neighboring segments with spin values I and m. 

Next, we consider the scattering probability of worms 
at a vertex with no kink on it. Suppose a spin- lowering 
worm hits the lower-left leg of the vertex from below in 
the coarse-grained picture. In order for this worm to 
be scattered, the worm and one of the legs of the vertex 
must be mapped onto the same line by the inverse map 
(Fig.^). There are two such cases: the case where the 
spin value is 1 [Fig.^(b)] on the legs on the other line and 
the case where it is [Fig.||(c)]. In the first case there 
are m different choices of the line, whereas we have m 
choices in the second case. Each individual choice in the 
first case has the weight pf^ ^ whereas that in the second 

(ss) 

case has the weight Piq ■ Therefore, the probability for 

choosing the first case is mp'fi^ / pim, whereas that for 

the second case is fhp'f^'^ / pim- If we choose the first case, 
the probability with which the worm is scattered in the 
direction F is 

P(=^^) {T\\- \ ). 

The worm scattering probability for the second case is 
given similarly. All in all, the probability of a spin- 
lowering worm being scattered into the direction specified 
by the directed graph F {^]) becomes 



P F 



m 
m 



, (ss) , _ (ss) J (ss) T_ (ss) 



.(3) 



The probability for going through (F —]) is simply equal 
to 1— (probabilities of the three proper scatterings). The 
symbol 



m 
m 



denotes the state where the spin states on the upper-left, 
upper-right, lower-left, and lower-right legs are l',m'J, 
and TO, respectively, and there is an incoming spin-raising 
(-I-) or spin-lowering (— ) worm on the lower-left leg. 
p(ss) (r|S) is the probabihty of the worm being scattered 
into the direction F when the initial state of the vertex is 
S in the split-spin representation. This probability coin- 
cides with that in the S = 1/2 case with the replacement 
Hp — > h. The scattering probability of a spin-raising 
worm can be obtained in the same fashion. 

The scattering probability at a vertex with a kink is 
simpler than that for a vertex with no kink on it, because 
in this case there is at most one type of vertex that may 
lead to a proper scattering (diagonal, horizontal, or turn- 
ing back). For example, suppose a particle jumps from 
left to right at the kink at the imaginary time r, and 
the spin-lowering worm is approaching the vertex on its 
lower-left leg. The local state S is given by 



S = 



I - 1 

l- 



m 
m 



1 



Then, the vertex's lower-left leg must be footed on pos- 
itive segments [^^^^(t) = 1] because otherwise no par- 
ticle can hop to the neighboring site there. Similarly, 
the lower-right, upper-left, and upper- right legs must be 
footed on negative, negative, and positive segments, re- 
spectively. There are Ifh such choices of segments, and 
all the choices are equally probable. Among them, there 
are m choices where the lower-left leg is footed on the 
segment where the worm is located. Therefore, the prob- 
ability of the worm being located on one of the legs of 
the vertex is m/{lfh) = . Then, the scattering prob- 
ability for F's corresponding to proper scatterings [i.e., 

r=T,/,i] is 



P F 



1 m 

m 



1 



1 

1- 



(4) 



for spin-lowering worms. Probabilities for spin-raising 
worms can be obtained similarly. 

Thus, we have described the way we derive the density 
of vertices and the scattering probability of worms from 
an algorithm in the split-spin representation. Although 
our description above may seem to give an actual proce- 
dure for coarse-graining mapping and its inverse, we do 
not perform these mappings in real simulation. They are 
only for deriving the density (^ and the probabilities (|^) 
and (^). In the actual simulation, we manipulate only 
coarse-grained variables. 

In order to complete the description of the algorithm, 
we have to specify the procedure for the pair creation 
and annihilation of worms. Again, this can be done by 
the coarse-graining map and its inverse. In the split-spin 
representation, the pair creation of worms is done sim- 
ply by choosing a point ((i,/j,),r) in the system with a 
uniform probability distribution. If there is a hole at the 
chosen point [i.e., ni^(T) = 0], we create a pair of spin- 
raising worms there. If there is a particle instead, we 



create spin-lowering worms. When coarse grained, this 
procedure is mapped to choosing a point from the whole 
space-time with uniform probability distribution and cre- 
ating a pair of spin-raising worms with probability 1/ {2S) 
or spin- lowering ones with probability 1/{2S), where / is 
the spin state at the chosen point. 

The moving worm travels according to the scattering 
process described above until it comes back to the origi- 
nal position ((z, /i), t) where the other worm waits. When 
coarse grained, this "coming-home" event is mapped to 
an event in which a worm comes back to How- 
ever, several other split-spin events are mapped to this 
same coarse-grained event. Namely, there are cases where 
the moving worm comes to the point corresponding to a 
different a spin, i.e., {{i,i>),T) with v ^ fx. Worms in 
this case should not annihilate. It has to be mapped, 
therefore, to a "going-through" event. Suppose that the 
worms are spin-lowering ones and that the local value 
of the coarse-grained spin is I (before the passage of the 
worm). Then, there are I cases in total which are mapped 
to the same coarse-grained state. Only one of them leads 
to the collision of two worms. Therefore, the probabil- 
ity of pair annihilation is and that for going through 
is 1 — Z""'^. For the same reason, the probability of an- 
nihilation should be if the worms are spin-raising 
ones. 

The whole procedure of one Monte Carlo sweep (MCS) 
with the algorithm described in this section can be sum- 
marized as follows. 

Step 1: Place vertices at random with the density, 

that depends on the local spin state, S. Set 

^count — 0- 

Step 2: Increase ricount by 1. Choose a point in the 
whole space-time at random, and create two worms 
there, one is to move and the other is to stay. For 
the moving worm, choose the initial direction of 
motion, upward or downward, with the probability 
1/2. Then choose its initial type, spin-lowering or 
spin-raising, with the probability 1/{2S) or 1/{2S), 
respectively. 

Step 3: Let the moving worm go until it hits a vertex or 
comes back to the original position where the other 
worm stays. If it hits a vertex before it comes back 
to the original position, go to step 4. Otherwise, go 
to step 5. 

Step 4: Choose the scattering direction F with the prob- 
ability P(F|I]), where E is the local spin state at 
the vertex before the worm's arrival. Change the 
type of the worm as specified by F. Then, go back 
to step 3. 

Step 5: If the moving worm is a spin-lowering one, let it 
go through the original point with the probability 
l — and go to step 3. Otherwise let it go through 
with the probability l — and go to step 3. If the 
moving worm does not go through, let it annihilate 
with its partner and go to step 6. 
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FIG. 4: Some examples of the final state E"" of scattering 
for which the initial state is E and the outgoing direction of 
the worm is specified by the directed graph F. A solid tri- 
angle denotes a spin-raising worm whereas a open triangle a 
spin-lowering one. The dots in the directed graphs represent 
segments for which spin variables are not changed by the scat- 
tering. The symbols in the parentheses are abbreviated forms 

of r. 



Step 6: If ricount is smaller than nmax, go back to step 2. 
Otherwise, erase all the vertices with no kink and 
go back to step 1. 

One Monte Carlo sweep is defined as a process between 
two successive resets of vertices (i.e., two successive pas- 
sages of step 1). The number n„iax, the number of pair 
creations of worms during 1 MCS, can be an arbitrary 
positive integer. We choose it so that every vertex may 
be visited by a worm once in average during 1 MCS. 



IV. THE XXZ MODELS 

Since the XXZ models are of particular importance, 
we summarize the probability of vertices and the scat- 
tering probability of worms for the models in Table |. 
Besides the coupling constants, the scattering probabil- 
ity depends upon the initial configuration of the scat- 
terer (i.e., vertex), the type of the worm ("spin-raising" 
or "spin-lowering"), and the incoming and outgoing di- 
rection. Because of the mirror image symmetries with 
respect to the horizontal and vertical axes, scattering 
probabilities for any two cases which can be transformed 
to each other by mirror image transformations should be 
the same. Therefore, without loss of generality, we as- 
sume that the incoming worm is located on the lower-left 
leg of the vertex. Then, the initial states can be catego- 
rized into six classes, each specified by the spin states on 
all the legs and the kind of the incoming worm (Table |). 

There are only a few possible final states for each initial 
state. Those final states can be specified by the outgoing 
direction (F) of the scattered worm (Fig.^). There are 
four such directions: turning-back, diagonal, horizontal. 
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FIG. 5: The six regions in the parameter space for the XX Z 
model with general S. The same as the one in Syljuasen and 
Sandvik's paper for S = 1/2. 



and straight, as indicated in the top row in Fig.^. The 
probabihtics for scattering in these directions are denoted 
by P(i |S),-P(^ andP(T |S) respectively, 

where E is the local state in the coarse-grained spin rep- 
resentation. In Table ^, we present the first three only. 
The probability for going straight P(t |S) can be readily 
obtained by subtracting the other three from unity. 

The scattering probability also depends upon the cou- 
pling constants, J and J'. From the algorithmic point 
of view, the whole parameter space is divided into six 
regions (Fig.||). Within each one of the six regions, the 
scattering probability is a simple analytic function of the 
coupling constant, and it is continuous in the whole pa- 
rameter space. However, its first derivatives are discon- 
tinuous at the boundary between two adjacent regions. 
In the case of S* = 1/2, the division is the same as that 
in the previous paper pO|. 



V. EFFICIENCY 

It is practically impossible to evaluate the efficiency of 
the algorithm for all possible combinations of coupling 
constants, the external field, the magnitude of spins, and 
the number of dimensions. Therefore, here we only show 
an example and make a few remarks concerning the effi- 
ciency of the algorithm described above. 

Of particular interest is the algorithm in region III, 
because the primary motivation for developing the al- 
gorithms based on the SS scheme is to solves the 
freezing problem of the conventional loop algorithms 
in this region. For S — 1/2, good performance was 
demonstrated ||T^ in the isotropic case \J'\ = J for var- 
ious values of H . Most importantly, no severe freezing 
was observed at low temperature. 

In what follows, we show that the present algorithm 
solve this problem for an arbitrary S. Several other di- 
rected loop algorithms (algorithms 1 - 4) in the origi- 
nal spin representation are also examined for comparison 



with the present algorithm. Algorithms 1-3 are ob- 
tained by tuning solutio ns o f the weight equation([A^) 
and the detailed balance(A4) so that the turning-back 
probabilities may be minimized. All of these three algo- 
rithms have exactly the same turning-back probabilities. 
Algorithm 1 is characterized by the vanishing probabil- 
ity for the diagonal scattering when the field is zero [i.e., 
lim/i^o -P(/' 1^) = for all S], whereas it is finite even 
at /i = in algorithm 2. Algorithm 3 is a mixture of 
algorithms 1 and 2. The details of these algorithms are 
given in Appendix Algorithm 4, on the other hand, 
is the heat-bath-type algorithm that can be obtained in 
the most straightforward way, although this is also a so- 
lution of Eqs. (AS) and (|A4|). In algorithm 4, one of 
the four possible directions F is chosen with the proba- 
bility proportional to the weight of the final state. To be 
specific. 



p(r|E) 



where E is the initial state of the vertex and YF is the 
final state of the vertex when the worm is scattered into 
the direction F. 

In order to check the validity of these algorithms, we 
first performed simulations for a small one-dimensional 
system [L — 4) and compared the results with the exact 
solution for various set of parameters, J', and (3. It 
turned out that all the algorithms yielded correct results 
with 1% or less of the statistical error. The present algo- 
rithm and algorithm 1 yielded roughly the same magni- 
tude of error whereas the other three yielded larger errors 
than the first two. 

For a longer chain {L — 64), 50 sets of simulations were 
performed using each algorithm where each set consists of 
20000 creations and annihilations of worm pairs. We can 
see the performances of five algorithms in Fig.^. Plotted 
in Fig.|is A{M^)Ny^^/L, where A(M^) is the estimated 
statistical error of the squared staggered magnetization, 
L is the system size and Ny is the average number of the 
vertices visited by the worm during its lifetime. Since 
the scattering process is the most time-consuming part 
of the code, the total CPU time is roughly proportional 
to the total number of scattering events of worms, in- 
cluding the "straight" scatterings. Therefore, the CPU 
time is proportional to Ny. This is why the statistical 

1/2 

error should be multiplied by Ny in order to make the 
comparison fair. In Fig.^, we can clearly see that the 
present algorithm performs as well as the best algorithm 
among the the other four (i.e., algorithm 1). Obviously, 
there is no exponential slowing down for the present al- 
gorithm and algorithm 1, as was the case with Syljuasen 
and Sandvik's algorithm for S = 1/2. 

For a larger system (L — 64) with zero magnetic field, 
the results of the present algorithm and algorithm 1 
agree with each other, while those of the other algorithms 
do not (Fig.^). We consider that the correct result for 
L = 64 is the one that is obtained by the present algo- 
rithm and algorithm 1, and that the other algorithms fail 
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TABLE I: The coarse-grained algorithm for the XX Z spin models. The density of vertices p, and the scattering probabilities 
of worms P. h = Hp/{2S), I = 2S - I, a,nd m = 2S - m. 



Region I Region II Region III 



Region IV 



Region V Region VI 



I rn 

I III 



A 



B 



B 



B 



A 



C 



fh{-J+J'-h) 
2C 
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2C 



I m 
1+ m 



P{1 |S) = 
Pi/- |S) 





m(J+j'-h) 
iA 

ra{J—j' — h) 
4/> 



i(-J-J'-h) 
2B 



mj 

2M 






i(J-j'-h) 
iB 



m(J+j'-h) 
4A 



I m 
l~ m 



Pi/- |S) 

Pj^ IS) 





■m(J+j'+h) 
iA 

m(J-j' + h) 
iA 



-i(-J-J'+h) 
2B 



mJ 



fh(-J-J'+h) 
2B 

m(J+j'+h) 
4S 

m J 

2B. 



m(-J + J' + h) 
+fh{-J - J' + h) 



mJ 
2B 



t(- J+j'+fc) m(- J+j'+fe) 
2A 2C 



mJ 
2A 

m(J-j'+h) 
iA 



mJ 
2C 



Z + 1 m 
/+ m + l 



P(.i IS) = 
P{/ |E) 

|s) 



J+j' + h 
1-2J 

J-j'-h 
LJLZ 





J + j' + h 
1-2J 

J-j'-h 
LJLZ 



Z — 1 m 
l~ m — 1 



Pil |S) = 
Pi/- |S) 

|S) 





J+j'-h 
1-2J 

J-j'+h 
L2J. 





J+j'-h 
1-2J 

J-j'+h 
L2J. 



Z + 1 m 
r m + l 



pa |E) = |E) = |E) = 0, and P(T |S) = 1 



Z — 1 m 
1+ m- 1 



pa |E) = P(/' |E) = PH |E) = 0, and Pa |S) = 1 



A = -[lm{J + J' + 3h) + {lfh + lrn){J-J' + h)+lfh{J + J'-h)] 

j/ I L -I 

B = lmh+{lm + Im) ^ , C = - [lm{J' + h) + lm{J' - h)] 



to achieve equilibrium within the performed simulation. 
We performed similar simulations for various values of H 
ranging from H = 0.0 to H — 2.0. It turned out that 
the good algorithms (the present and algorithm 1) always 
perform better than the bad ones (algorithm 2-4) al- 
though the difference between them becomes smaller as 
the field is increased. 

We can explain these facts in terms of the compati- 
bility of the clusters with the order parameter. One of 
the reasons why the conventional loop algorithm works 
well even in the vicinity of the critical point is the ac- 
cordance between the typical cluster size and the cor- 
relation length. This accordance is guaranteed by the 
two features of the algorithm: (i) independent flipping 
of clusters, and (ii) perfect ordering within each clus- 
ter. Although the staggered magnetization is not strictly 
the order parameter in one dimension, this criterion of 
good performance of loop-cluster algorithms still applies 
because finite but relatively long-range correlation ex- 
ists even in one dimension. It is easy to sec that the 
present algorithm and algorithm 1 satisfy both condi- 
tions (i) and (ii) in the zero-field limit whereas the other 
algorithms do not satisfy condition (ii) regardless of the 
field. This is the reason why the former two algorithms 



perform much better than the latter three in the weak 
field region. Therefore, we expect that the difference in 
the efficiency is even more pronounced near real phase 
transitions such as these in three-dimensional systems. 

The results of the five algorithms, all based on the SS 

scheme, illustrate that it is not trivial in general to ob- 
tain the most efficient algorithm among many possible 
ones and also that the straightforward heat-bath algo- 
rithm is rather poor in some important cases. It should 
be noted that the coarse-grained algorithm discussed in 
the present paper satisfies criterion (ii) mentioned above 
for an arbitrary XXZ model when the external field is 
vanishing. Therefore, we consider that the present algo- 
rithm performs well for a rather wide class of models in 
the weak magnetic field region. Even for strong magnetic 
field, we consider that the present algorithm is at least 
as good as most of the other algorithms based on the SS 
scheme as we see above. 

For the performance of the present algorithm in the 
regions other than III, we cannot conclude much at the 
moment. In the case of >S = 1/2, the algorithm in region 
I is equivalent to the Wolff version (i.e., the single-loop 
update version) of the loop algorithm when the exter- 
nal field is vanishing. A good performance of the al- 
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in the classical Ising limit) the conventional loop-cluster 
algorithm works efficiently. For region II, it is not known 
whether this is a real problem or not. The algorithm 
in region IV is reduced to a single-spin-flip Metropolis 
algorithm in the limit of J'/ J — > and h ^ oo. The 
performance of the Metropolis algorithm should be good 
in this limit, although the region is physically not very 
interesting. 



VI. SUMMARY 



FIG. 6: The statistical error in the estimate of the squared 
staggered magnetization multiplied by the square root of the 
average number of scattering events during the lifetime of a 
worm. The system is the S = 1 antiferromagnetic Heisenberg 
chain of length L = 64 with a uniform magnetic field H — 0.1. 
Each point is a result of 50 sets of simulations, where each set 
consists of 20000 pair creations and annihilations of worms. 
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FIG. 7: The squared staggered magnetization estimated with 
the present algorithm and those with the other algorithms for 
the S — 1 antiferromagnetic Heisenberg chain with L = 64 
and H — 0. The rest of the condition for the simulation is 
the same as that for the previous figure. 



gorithm in this region has been demonstrated in many 
applications 0, as well as for 5 > 1/2 ||, Since the 
turning-back probability is vanishing in the whole region 
I, we expect good performance of the algorithm not only 
on the line h — but also in the whole region I. The al- 
gorithm in regions II and VI, on the other hand, may not 
be very efficient because of the presence of the relatively 
large turning-back probability. It is easy to see that, in 
the classical limit {J'/J — *■ oo) at zero external field, the 
turning-back probability dominates in regions II and VI, 
leading to poor performance. For region VI, this may 
not be very problematic because in this region (at least 



We have proposed an algorithms based on the 
Syljuasen and Sandvik scheme by introducing the split- 
spin representation and the coarse-graining procedure. 
The algorithm is a natural extension of the directed loop 
algorithm for S — 1/2, in that the present algorithm co- 
incides with it for S = 1/2. In addition, the present 
algorithm is a natural extension of the conventional loop 
algorithm, because if the external field is vanishing, the 
present algorithm can be obtained through coarse grain- 
ing the conventional split-spin loop algorithm. 

Compared with the algorithms in the split-spin rep- 
resentation, the coarse-grained algorithm requires much 
smaller memory, in general. In particular, when the al- 
gorithm consists of vertices with more than four legs, as 
is generally the case with the loop algorithms for mod- 
els with high order interaction terms, the memory can 
be reduced considerably. The coarse-grained algorithm 
is also advantageous since the codes based on it can be 
very easily modified for other models (we only need to 
change the arrays of the probability tables). 

In the case of the S — 1 Heisenberg antiferromagnet 
in one dimension, the algorithm's performance is almost 
the same as the best algorithm obtained by directly work- 
ing on the original spin representation. Many other al- 
gorithms can also be obtained in the same way. How- 
ever, most of them, including the heat-bath algorithm, 
are much worse than the present one. Ex iste nce of algo- 
rithm 1, a good direct solution to Eqs. (A2) and (K4), 
suggests the existence of similar solutions for an arbitrary 
S. We have not succeeded in finding a complete set of 
such solutions, although we believe that such solutions 
exist. It would be an interesting future problem to find 
such solutions for various models. 

The coarse-graining procedure presented in this paper 
applies not only to the XXZ spin systems but also to 
any model for which a directed loop algorithm can be 
constructed in the split-spin representation. For exam- 
ple, on-site easy axis or easy plane anisotropy terms may 
be treated as the couplings between a spins on the same 
site. Another example is the SU(N) model where the 
split-spin algorithm is known. 
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APPENDIX A: GENERAL FORMULATION OF 
THE SS SCHEME 

In general, a directed loop algorithm is characterized 
by the density of vertices /9(S) and the scattering prob- 
ability of worms P(r|S). The density p{l,m) is simply 
given by 



I m 

ni 



(Al) 



where the weight M^(S) is defined as 



W 



V m' 



{ch'Smm' - {l',m'\Hij\l,m)) 

xA(0 < / + e < 25), (A2) 



where A(" • ■ • ") = 1 when " • ■ • " is true and otherwise. 
The symbol e stands for the integer by which the worm 

changes the spin value, e.g., e — —2 for a ( ) worm 

that lowers the spin value by 2. The variable c is the 
only free parameter related to the vertex density. 

The scattering probability, on the other hand, has a 
lot of freedom. The algorithm can be explained very 
clearly by introducing an extended weight VF(E,r) that 
is related to the weight W^(I]) as 



w^(i])-5]w^(s,r). 



(A3) 



Here, we consider a scattering event in which the initial 
state of the vertex is S and the out going direction of 
the worm is F. We denote the final state of this event as 
Y,^ . It should be noted that the 'state' S is directed in 
contrast to the state in the ordinary Monte Carlo sim- 
ulation. We consider the balance between an arbitrary 
sequence of scatterings and its reverse. Each sequence 
starts from pair creation of the worms and ends at pair 
annihilation. The detailed balance condition should be 
considered between S and the reverse ofE^. Therefore, 
the detailed balance condition is expressed as 



VF(i],r) = VF(sr,r) 



(A4) 



Here, E'" is the reverse of E'" obtained by inverting the di- 
rection and changing the type of the worm in Y.^ , whereas 
f is the inverse of F, obtained by inverting the direction 
of the arrow in F. 

It is worth mentioning that Eqs. (A_3) and (A4) are 
quite similar to the equations that appear in the general 
formulations of the loop-cluster algorithm B Oj. The 



only difference is that here we use directed graphs and di- 
rected states whereas only nondirected graphs and states 
appear in the conventional loop-cluster algorithm. It is 
easy to see that in the case of zero magnetic field the 
extended weight in the present scheme can be made in- 
dependent of the directions of states and graphs, and all 
the equations discussed in this appendix coincide with 
those for conventional loop-cluster algorithms. 

Once we obtain any set of constant c and po sitive 
VF(I],F)'s that satisfy Eqs. (|a|), (^, and (^), we 
obtain a scattering probability P(F|S) of worms that sat- 
isfies the detailed balance condition as 



P(F|I]) 



(A5) 



APPENDIX B: SOME DIRECTED LOOP 
ALGORITHMS FOR THE S = 1 
ANTIFERROMAGNETIC HEISENBERG MODEL 



In Tables || and |T|, we show the weight and the ex- 
tended weight that satisfy Eqs. (|A^), ([A3), and ( [A^ ) for 
the S = 1 antiferromagnetic Heisenberg model. The ver- 
tex density and the wo rm-s catter ing probability can be 
obtained through Eqs. ( |Al| ) and ( A5). 

Table || shows the scattering probability for worms 
that change values of spin by 1. The weights contain two 
free parameters A = J — B and A' — J — B' . Note, how- 
ever, that the turning-back probability does not depend 
on the choice of the free parameters. Algorithm 1-3 
correspond to the following choices, respectively: 

algorithm 1, A= A' ^ J, 

algorithm 2, A = J, A' = 

algorithm 3, A = 0.9J + 0.1^, A' = 0.9^ +0.1J. 



In addition to the worms changing spin values by 1 , we 
can introduce worms that change spin values by an arbi- 
trary amount. For S = 1, we can introduce ±2 worms. 
In the examples presented in Sec. ^ we used both ±1 
worms and ±2 ones. When a pair of worms are created, 
the type of worm is chosen with equal probability from 
all possible ones. For instance, when the initial spin state 
is Z = 2 at the point chosen for the pair creation, a — 1 
or —2 worm is possible. Each one of them is chosen with 
probability 1/2. The extended weight for ±2 worms are 
listed in Tabic 0. 



In fact, only the ±1 worms are necessary for making 
the algorithm ergodic and for satisfying the detailed bal- 
ance. Although it is likely that the ±2 worms are use- 
ful for improving the efficiency of the algorithm for more 
complicated models, we have not encountered such a case 
yet. 
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TABLE II: A ±1 worm solution (VF(E,r)) to the detailed balance equation for the 5 = 1 antiferromagnetic Heisenberg model. 
Applies only if < Hp < 4J. Free parameters A, B, A' , and B' are related to each other hy A + B = A' + B' = J. 
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TABLE III: A ±2 worm solution to the detailed balance equation for the S = 1 antiferromagnetic Heisenberg model. Applies 

oulv if <- < 2.J. 
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